Quantum computing and preconditioners for hydrological linear systems

Modeling hydrological fracture networks is a hallmark challenge in computational earth sciences. Accurately predicting critical features of fracture systems, e.g. percolation, can require solving large linear systems far beyond current or future high performance capabilities. Quantum computers can theoretically bypass the memory and speed constraints faced by classical approaches, however several technical issues must first be addressed. Chief amongst these difficulties is that such systems are often ill-conditioned, i.e. small changes in the system can produce large changes in the solution, which can slow down the performance of linear solving algorithms. We test several existing quantum techniques to improve the condition number, but find they are insufficient. We then introduce the inverse Laplacian preconditioner, which improves the scaling of the condition number of the system from O(N) to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\sqrt{N})$$\end{document}O(N) and admits a quantum implementation. These results are a critical first step in developing a quantum solver for fracture systems, both advancing the state of hydrological modeling and providing a novel real-world application for quantum linear systems algorithms.

Quantum algorithms for solving linear systems offer a potential for a practical quantum advantage over classical algorithms 1 . However, several technical and conceptual challenges remain before a meaningful, real-world example of quantum advantage can be exhibited for linear systems. On the technical side, many of these quantum linear systems (QLS) algorithms can only be implemented on future error-corrected hardware 1 , while others are designed for today's noisy intermediate-sized quantum devices but feature a less pronounced speed-up 2 . On the conceptual side, existing algorithms feature numerous caveats that must be satisfied in order for good performance 3 . The ideal system Ax = b for exhibiting quantum speed-up satisfies: (1) A is well-conditioned, (2) the quantum operator e −iAt can be efficiently implemented, (3) a quantum state proportional to b can be efficiently prepared, (4) only specific values or statistics of the solution x are of interest, (5) there does not exist a classical algorithm which exploits these or other properties of the system to solve it equally as fast. The final point is particularly notable, as large systems of equations can often be replaced by small systems of equations which retain sufficient accuracy (and can be solved classically with modest cost) 4 . Existing applications for QLS algorithms have therefore been largely synthetic or specific examples, carefully chosen to avoid these constraints but also lacking in real-world application 5,6 .
In this work we initiate the study of simulating fluid flow through fracture systems with QLS algorithms. Hydrological flow is a very challenging problem in geophysics, because the contrast between the scales on which simulations are often done (kilometers or larger) and the scale of heterogeneities (centimeters or smaller) requires discretizing the problem over very large meshes. Here we address the common problem of determining the pressure of a subsurface liquid (e.g. water or oil), either at a specific location, i.e. at a well, or averaged across a wide region, i.e. effective permeability. Extracting the pressure at an individual location or studying averaged properties (rather than needing exact pressure at every point in the system) addresses point 4 in our list of QLS requirements. As we argue in section "Fracture networks and coarse-graining", this type of problem cannot be reduced to a smaller system of equations without losing critical information, and so current classical techniques cannot capture the full scope of the problem at real-world scale (point 5). Since the matrices A appearing in subsurface flow applications are generally sparse, and the vectors b are relatively uniform, points 2 and 3 are plausibly addressed via quantum RAM 3 . However, these points are quite nuanced and will be studied in detail in a future work.
Our focus for this work is point 1, improving the condition number of the linear system. The condition number of a matrix describes how sensitive the system is to small perturbations, or more specifically, how much an error in the right hand side, b, propagates to an error in the solution, x, in the worst case. Easy matrices, such as www.nature.com/scientificreports/ unitary matrices, have condition number 1 and increasing the condition number makes the linear system harder to solve. Both quantum and classical linear systems algorithms generally perform worse on systems with large condition numbers. Classical methods to reduce the condition number, known as preconditioners, have been successfully developed and tailored specifically for fracture systems, the most prominent example being multigrid preconditioners 7 . However, they cannot simply be ported directly onto quantum computers due to the drastically different underlying technologies and algorithmic design constraints 8 . Therefore, in order to use quantum computers to solve large linear systems for fracture networks of relevance and outpace classical techniques, an effective quantum preconditioning algorithm for fracture networks must be developed. A good preconditioner is generally tied to the type of linear system being solved and exploits some underlying mathematical or physical structure in the problem. Meanwhile, one must also find an efficient way of calculating the preconditioned form of the system, which will be necessary to do on the quantum computer itself for very large systems due to memory constraints. These two constraints-tailored to a specific application while also having an efficient quantum implementation-are likely to play a prominent role in future attempts to solve large, interesting linear systems on quantum computers. We show that existing techniques to reduce the condition number either do not have quantum implementations or are not well-suited for fracture systems. However, we find that the inverse Laplacian is both effective on fracture systems and can be efficiently implemented on a quantum computer. See Fig. 1 for a summary of our results. Furthermore, we show that realistic fracture systems can be solved by a quantum computer with a polynomial advantage over classical approaches. Finally, we discuss the remaining algorithmic bottlenecks which need to be resolved to unlock the full potential of QLS for fracture systems.

Results
In this section, we describe the geophysical problem in further detail and argue why it is not fully tractable for classical computers (section "Fracture networks and coarse-graining"). We then review the three existing quantum preconditioner algorithms and study how well they improve the condition number of several synthetic fracture network examples. The condition number (without preconditioning) of these examples all scale as O(N), and the best of the preconditioners reduces this scaling to roughly O( √ N) (section "Test of existing quantum preconditioners"), which is potentially sufficient to produce a polynomial advantage over classical techniques. Finally, we analyze how incorporating the preconditioner into a full QLS algorithm might affect overall runtime (section "Algorithmic scaling with inverse Laplacian preconditioner").
Fracture networks and coarse-graining. In a complex fracture network, fractures of many scalesfrom kilometers to centimeters-intersect. Critically, small fractures cannot generally be neglected because these can transform the network topology radically, e.g. pushing a system over a percolation threshold, see Fig. 2. Small fractures may also collectively contribute a large surface area to the network providing a critical connection between fractures and the underlying rock. When modeling flow in these networks, it is therefore critical to include the full range of fracture scales, which has led to the development of advanced meshing www.nature.com/scientificreports/ techniques 9 and high-performance simulators 10 . However, these approaches do not provide a viable path to modeling the full range of scales. Even state-of-the-art high-performance computers and cutting-edge methods can only model large fracture networks with fracture lengths varying over three orders of magnitude 11 . Modeling real-world fracture networks to a high degree of accuracy requires meshes far beyond current or future classical capabilities-e.g., a 1km domain with a 1cm resolution would require 10 15 degrees of freedom. Larger systems of equations would be needed for larger domains or more finely resolved meshes. Numerical models of subsurface flow in a fracture network are based on a discretized version of 7 where k is the permeability, f is the fluid flux, and h is the pressure. Each of k, f, and h is a funciton that varies spatially throughout the domain. The pressure is key for applications of subsurface flow such as waste fluid disposal and hydraulic fracturing 12 . It is also critical for transport applications, since pressure gradients drive the transport. In fractured systems, k is highly heterogeneous with a sharp increase when moving from the rock matrix (where k is small) to a fracture (where k is large). We discretize Eq. (1) using a two-point flux finite volume method, which is one of the standard numerical schemes used in subsurface flow solver code bases including: fehm 13 , tough2 14 , and pflotran 15 . The two-point flux finite volume method ensures mass conservation, which is a highly desireable property for these numerical solvers. This results in a coefficient matrix which is sparse, symmetric, and positive definite. We treat the permeability of the rock matrix as being constant. While the approximation is imperfect, it is a major step up from discrete fracture networking approaches which effectively treat the rock matrix as having zero permeability 7 . This approach does, however, capture the basic physics of fracture flow-most but not all of the flow occurs in the high-permeability fractures.
In this work, we study a variety of 2D fracture network models. The simplest system we studied involved two fractures intersecting in a -configuration, and we then studied fractal-style recursion of the -system to generate more complicated fracture networks, see Fig. 3. The relative permeability of the fractures as compared to the underlying rock is a critical parameter in the analysis of fracture systems, and we studied five different types: Fracture Type 1: "Simple, Low" is the -system with fracture 10% more permeable as the underlying rock. Fracture Type 2: "Simple, High" is the -system with fractures 10 4 times more permeable than the underlying rock. Fracture Types 3 and 4: "Fractal, Low/High" are the same as above, but with the fractal system. Fracture Type 5: "Fractal, Var. " is the fractal system where the fractures have permeability contrast that scales down as the fractures get smaller, i.e. largest fractures have contrast 10 4 and the smallest fractures have contrast 1.1, with the contrast scaling down as the (fracture length) 1/2 , commonly used in pracice 11 .

Test of existing quantum preconditioners. Recently developed QLS algorithms provide a novel path
to modeling the full complexity of fracture networks. Since a quantum computer with n qubits can represent a 2 n dimensional vector, vast systems of equations can be solved with a small number of qubits. That is, for the 1km domain with a 1cm resolution problem described above, a quantum computer could require as few as O(log(10 15 ) ≈ 50) qubits, whereas a classical computer would require O(10 15 ) classical bits. While additional www.nature.com/scientificreports/ quantum resources will likely be necessary, either in the form of ancilla qubits or quantum RAM, this highlights the raw potential of the quantum approach for the huge systems of equations necessary to accurately model realistic subsurface flow. Furthermore, the computational complexity of quantum linear systems algorithms can in some cases be exponentially better than the best classical counterparts 16 .
The original QLS algorithm introduced by Harrow, Hassidim, and Lloyd (HHL) 1 solves a sparse N-variable system of equations Ax = b with a runtime of O(log(N)κ(A) 2 ) , where κ(A) := �A��A −1 � is the condition number of A (in this work we use A with no subscripts to refer to the 2-norm (or operator norm), i.e. �A� ≡ �A� 2 , and explicitly use A F to refer to the Frobenius norm). The best classical algorithm, the Conjugate Gradient method, runs in O(N √ κ(A)) on sparse matrices, so the quantum algorithm provides an exponential speed-up when κ(A) is small (for notational clarity, we often use ' κ ' alone to signify κ(A)).
Since the introduction of the HHL algorithm, many more QLS algorithms have been introduced, including improvements to HHL 16,17 , adiabatic approaches 18,19 , and variational algorithms implementable on near-term quantum computers 2,5 . A feature common to all of these quantum approaches is a scaling with κ that is linear at best, as compared to the O( √ κ) scaling of the best classical approaches. A common technique in classical analysis is to further reduce κ by using a preconditioner. This technique relies on finding a matrix M such that κ(MA) ≪ κ(A) . One then finds the x satisfying MAx = Mb . The matrix M is generally dependent on the specific matrix A, and different preconditioning approaches have been developed for different contexts.
Despite the significant interest and activity in QLS algorithms, relatively little work has been done to develop application-specific preconditioning algorithms. It is important to emphasize that in this work we only consider preconditioning algorithms which are implementable on a quantum computer, as opposed to any sort of hybrid classical-quantum preconditioning method. While we do not rule out the possibility of such an approach, the extreme memory requirements of full-scale subsurface flow problems (as described above) suggests that calculating the preconditioned matrix classically would be intractable. There are currently three general purpose quantum preconditioning algorithms in the literature: the circulant method 20 , the sparse approximate inverse method 21 , and the fast-inverse method 22 . These algorithms are described in detail in section "Methods", here we give only the salient points.
The circulant method is a one-size-fits-all approach, that is, the only input is a matrix A and the output is a preconditioner M. With SPAI one gives a sparsity pattern for the preconditioner M, and several techniques have been developed for determining good sparsity patterns for fracture systems 23 . The fast-inverse method is designed for systems of the form A + B , where A −1 can be easily calculated. One then uses A −1 as the preconditioner. As discussed in section "Methods", fracture systems can be decomposed into + A F , where the Laplacian describes the system in the absence of fractures, and A F is the contribution of the fractures. Because the singular value decomposition of the Laplacian is known 24 , −1 can be efficiently calculated and used as a preconditioner.
In Fig. 4 we numerically evaluate the effect of the circulant, SPAI, and inverse Laplacian preconditioners on all of the fracture types described in section "Fracture networks and coarse-graining" up to O(10 6 ) variables. To estimate how the preconditioner performance scales with N, we perform a linear regression on the logarithm of the final four data points for each preconditioner applied to each fracture type. We do not include the first points as they sometimes exhibited variance due to small matrix sizes, however for N > 10 3 , all of the results show clear exponential scaling in the number of qubits and polynomial scaling in the number of equations. We find that the circulant and SPAI methods are poor choices for the fracture systems under consideration. While these two preconditioners consistently reduce the condition number of the system, they do not improve how the condition number scales in N, which is necessary to unlock the full potential of QLS algorithms for large problems. The inverse Laplacian preconditioner, however, does meaningfully improve the scaling of the condition number. In the cases with low permeability contrast, the condition number of the system −1 A is very low, scaling as ≤ O(N 0.05 ) . The high permeability contrast systems do not fair as well, with the κ of the preconditioned fractal system scaling as O(N 0.6 ) . The fractal system with variable permeability, which is the most realistic of the systems under consideration, has a preconditioned κ which asymptotically scales as approximately O(N 0.55 ). www.nature.com/scientificreports/ Algorithmic scaling with inverse Laplacian preconditioner. Identifying a preconditioner M which reduces the condition number of a system A is generally not sufficient to guarantee good performance of a QLS algorithm to solve the preconditioned system MAx = Mb . This is because one must find an efficient way to calculate the product MA in such a way as to make it readily accessible on a quantum computer. This can be accomplished either through a classical oracle which efficiently calculates MA, or through a quantum algorithm which uses oracles for A and M and then calculates MA. As described previously, in this work we only study cases where the preconditioner can be applied on the quantum computer itself, due to the considerable memory demands of full-scale hydrological simulations. Unfortunately, simply calculating MA in the general case on a quantum computer with generic matrix multiplication adds an O(N 2 ) overhead 25 and would remove any benefit from the reduced condition number. Each of the three methods previously discussed get around this limitation through clever techniques: the circulant method calculates the product with the quantum Fourier Transform algorithm, the SPAI method exploits sparseness of M and A, and the fast-inverse method assumes efficient block-encodings of M and A. As described in section "Test of existing quantum preconditioners", the data presented in Fig. 4 shows that the circulant and SPAI preconditioners do not meaningfully improve the scaling in N of the condition number for the fracture examples considered in this paper. Therefore, it is not worth determining the full algorithmic runtime for implementing either of these preconditioners in a full QLS algorithm. However, the inverse Laplacian reduces the condition number to scale between O(N 0.05 ) to O(N 0.6 ) , which is potentially sufficient for an advantage over classical algorithms. It is thus instructive to estimate the impact the inverse Laplacian preconditioner has on the overall runtime for solving subsurface flow systems.
While implementing the inverse Laplacian preconditioner into the original HHL algorithm is potentially possible, the QLS algorithm of Tong et al. 22 already gives a direct method of applying the preconditioner and solving the resulting system. We can therefore use the scaling of their algorithm as a proof-of-concept to gain an understanding of whether the inverse Laplacian improves the condition number sufficiently to recover some quantum speed-up. This analysis purposefully ignores intricacies resulting from points 2 and 3 in our list of QLS caveats, i.e. efficiently turning the classical data A and b into appropriate quantum states and operators. We emphasize that this analysis is intended as a conservative estimate of how a preconditioned QLS algorithm might scale when solving fracture systems. There are hopefully more efficient implementations, which we discuss more in section "Discussion" and will explore further in future work.
As we show in section "Methods", the fast-inverse QLS algorithm with −1 as the preconditioner gives a runtime bounded below by www.nature.com/scientificreports/ The scaling in N of �A − � and �A −1 � are dependent on the exact fracture systems being studied and must be determined experimentally. In Fig. 5 we show the scaling of these components for the fractal system with variable contrast. We focus on this particular example since it is the most realistic of the different fracture types. As was the case in section "Test of existing quantum preconditioners", we estimate the large-N scaling of the different components by linear regression (on the log-log plot) of the data points for N > 10 3 . In Table 1 we summarize the scaling of each component as well as the overall scaling (modulo log(N) ) compared with the scaling for Conjugate Gradient on the same systems.
Using the fast inverse QLS algorithm with the inverse Laplacian preconditioner, we can potentially achieve a polynomial improvement over the best generic classical scaling for all fracture systems considered here. This approach utilizes block encodings of −1 and A − to calculate the product −1 A . However, since � −1 � scales linearly in N, the block encoding takes at least this long. Future algorithms for QLS, specifically tailored to fracture systems, could be developed to calculate −1 A even more efficiently by exploiting the sparseness of A.

Discussion
In this work, we have initiated the study of using QLS algorithms to solve systems of equations describing subsurface flow. As we argue in section "Fracture networks and coarse-graining", capturing the full behavior of these systems at real-world scale is prohibitively memory-intensive for classical computers. Quantum computing provides an alternative path, potentially using significantly less resources and offering improved scaling in problem size. However, several conceptual problems must be addressed (in addition to the need for improved quantum computing hardware). The most prominent of these issues are the poor condition number of these systems and the means of loading information onto the quantum computer. Here we have studied the first problem through the use of quantum preconditioners, and we leave the latter problem to future work.
We have shown that two previously introduced quantum preconditioners, the circulant and SPAI methods, do not improve the scaling in N of κ(A) and therefore will not help gain a quantum advantage for these fracture systems. However, the inverse Laplacian is an effective preconditioner for fracture systems and readily admits a quantum implementation. In particular, it can be implemented via the fast-inverse QLS algorithm, and the overall scaling of this solver scales better than the best generic classical algorithm.   Table 1. Scaling of the various components entering the overall scaling of the fast-inverse QLS using −1 as the preconditioner. The Overall scaling reported is modulo log(N). www.nature.com/scientificreports/ In comparing against classical techniques, we have so far not addressed the fact that for PDE-based systems on a uniform mesh, such as those considered here, more specialized methods can be used. Geometric multigrid methods, which exploit the structured mesh, can solve systems of equations in O(N) or O(N log N) 26 . Due to the extensive caveats and stringent requirements attached to QLS algorithms, it is noteworthy (though by no means sufficient) that our preliminary results suggest that quantum performance will be at least comparable to stateof-the-art, highly tuned classical techniques. This is in addition to the memory requirements which classical techniques inevitably hit. Still, it is clearly necessary to further refine QLS algorithms with an even greater eye towards the specific physics and mathematical structures of the fracture systems at hand.

Scaling in N of:
An obvious area of improvement is a more efficient quantum means of implementing the preconditioned system −1 A . In the most realistic case we studied, the fractal system with variable contrast, κ(� −1 A) asymptotes to roughly O(N 0.55 ) . Therefore in principle the scaling of solving just the preconditioned system, e.g. with recent adiabatic QLS algorithms 18 , would be O (N 0.55 log N) , a significant improvement over the geometric multigrid methods.
Alternatively, while the inverse Laplacian opens the door for polynomial speed-up over classical, an even better preconditioner is needed for exponential speed-up. For example, while a direct quantum port of the classical multigrid methods is not plausible, some of the ideas may be used to construct an analagous approach that can be implemented on a quantum computer. Reducing the condition number scaling to O(log N) in the quantum context would then be possible.
This study shows that fracture networks are a challenging real-world problem with a potential for serious advancements from quantum computation. The large linear systems necessary to accurately model flow behavior are sparse yet cannot be coarse-grained. The condition number κ of the systems tends to scale linearly with N, however the inverse Laplacian preconditioner, which readily admits a quantum implementation, can improve this scaling considerably, and it is likely evern further advancements can be made in κ . Future work will be devoted to incorporating our application-specific preconditioning techniques into a full quantum linear solver, ideally targeting an implementation on NISQ devices.

Methods
In this section we provide a more detailed review of the preconditioner methods evaluated in Section "Test of existing quantum preconditioners", as well as a derivation of the scaling bound for the fast-inverse algorithm found in section "Algorithmic scaling with inverse Laplacian preconditioner". Circulant preconditioner. The circulant preconditioner method of Shao et al. 20 gives an efficient quantum implementation of a circulant preconditioner C based on the quantum Fourier transform F. An n × n matrix C is circulant if C ij = C (i−j) mod n . The use of circulant preconditioners in classical applications is motivated by the fact that, for a given circulant matrix C and an arbitrary matrix A, CA and C −1 A can be computed in O(n log n) steps using the fast Fourier transform. Circulant preconditioners are particularly useful in solving Toeplitz systems 27 .
For an arbitrary matrix A, one can construct the circulant preconditioner via where F jk = 1 √ n ω jk with ω = e −2πi/n . C −1 (A) is then used the preconditioner. F can be efficiently implemented via the quantum Fourier transform, and the middle term simplifies to An algorithm for efficiently preparing the state in Eq. 4 is given in 20 . This approach works for arbitrary dense non-Hermitian matrices, however there is no upper bound on κ(CA) , and in practice for random dense matrices κ(CA) = O(κ(A)).
Sparse approximate inverse. The Sparse Approximate Inverse (SPAI) approach for solving a system Ax = b attempts to find a matrix M such that MA ≈ I , where M has a (user-defined) sparsity pattern. For example, if one gives a sparsity pattern involving n non-zero rows and d non-zero elements per row, then M is given by solving n × d independent least squares problems. The trick with this approach is determining which sparsity pattern to choose for M.
Clader et al. 21 showed that the preconditioned system can be solved via a slightly modified version of the HHL algorithm. The overall scaling for actually solving Eq. (5) with error ǫ is In section "Test of existing quantum preconditioners" we adopt the relatively standard approach of using the sparsity pattern of A for M. One can also try other methods 23 , which can significantly reduce the condition number, but again do not improve the scaling in N for the fracture systems studied here. Finally, d does not need to be a constant. As long as κ(MA) = O(1) , the sparsity pattern can scale with d ∝ N ≤1/7 in order to at least recover  www.nature.com/scientificreports/ some quantum advantage. However, for the systems and sparsity patterns considered here, a small increase in d has a corresponding small decrease in κ(MA) . For example, when applying the technique described in 23 to the "simple" systems, i.e. fractal depth 1, increasing the density of M by a factor of five only decreases κ(MA) by a factor of two. It is therefore difficult to imagine a system size or sparsity pattern where such an incremental increase in d could produce sufficiently large reductions in κ(MA) as to make the procedure worthwhile.
Fast inverse. The fast inverse method of Tong et al. 22 solves linear systems Ax = b where A can be decomposed as where �A big � ≫ �A small � . They then give a QLS algorithm which uses A −1 big as a preconditioner and solves the system (I + A −1 big A small )x = A −1 big b with scaling bounded by �A small �, �A −1 big �, and �A −1 �. This technique is dependent on efficient block-encodings of A −1 big and A small . An (α, m, ǫ)-block-encoding of the matrix A is given by the unitary U A : where * denotes arbitrary matrix blocks, α is a rescaling constant such that �U A � = 1 , and the error ǫ is bounded by �A − α(�0 m | ⊗ I n )U A (|O m � ⊗ I n )� ≤ ǫ . Since the magnitude of α plays a critical role in the scaling of this algorithm, we note that �U A � = 1 implies that α ≥ �A�.
In order to use A −1 big as the preconditioner, A big must be fast-invertible. A matrix M is fast-invertible matrix if one can efficiently prepare a (�(1), m, ǫ)-block-encoding U ′ M of M −1 . This requires access to an oracle for M −1 , and the number of queries to this oracle in preparing U ′ M must be independent of κ(M) . For example, if M is normal, and the eigenvalue decomposition M = VDV † gives a V that can be efficiently implemented in a quantum circuit and the elements of D can be accessed through an oracle, then M is fast-invertible.
The fast-inverse QLS algorithm takes as inputs an (α s , m s , 0)-block-encoding U s of A small , and an They then use a modified version of the quantum singular value transformation 28 to construct a block encoding of (A big + A small ) −1 with error ǫ in applications of U s , U ′ b along with their inverses, controlled versions, and other primitive gates. Here σ min is a lower bound for the smallest singular value of I + A −1 big A small , i.e. the preconditioned system. This approach has the benefit of providing an upper bound on the condition number of the preconditioned matrix, with the downside of needing a decomposition of A that matches a lengthy list of requirements. For fracture problems, we have the natural decomposition of where the Laplacian describes the flow in the absence of fractures, and A F denotes the fracture matrix. Fortunately the discretized Laplacian is normal and has a known eigenvalue decomposition 24 , therefore it meets the fast-invertible conditions and we may use it as the preconditioner. However, we have no guarantee that � � ≫ �A − � , which is required to get good scaling. Still, we can numerically test the scaling of the algorithm to see how it performs in the absence of performance guarantees.
The parameters contributing to the performance of this algorithm, Eq. (9), are the block-encoding parameters α ′ b and α s , along with a lower bound on the smallest singular value of the preconditioned system MA, σ min . In order to assess the potential usefulness of this algorithm for our application, we will explore the most optimistic values for these parameters. Due to minor technical details, we rescale the entire system by � −1 � , which gives A big = � −1 � and α ′ b (the block-encoding parameter for A −1 big ) = �(1) . We also have A small = (A − �)�� −1 � , so α s ≥ �A − ���� −1 � , and 1/σ min ≥ �A −1 �� . Therefore the overall scaling for the fast-inverse QLS algorithm is bounded below by

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.